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The family of generalized-harmonic gauge conditions, which is currently used in Numerical Rela- 
tivity for its singularity-avoidant behavior, is analyzed by looking for pathologies of the correspond- 
ing spacetime foliation. The appearance of genuine shocks, arising from the crossing of characteristic 
lines, is completely discarded. Runaway solutions, meaning that the lapse function can grow without 
bound at an accelerated rate, are instead predicted. Black Hole simulations are presented, showing 
spurious oscillations due to the well known slice stretching phenomenon. These oscillations are made 
to disappear by switching the numerical algorithm to a high-resolution shock-capturing one, of the 
kind currently used in Computational Fluid Dynamics. Even with these shock-capturing algorithms, 
runaway solutions are seen to appear and the resulting lapse blow-up is causing the simulations to 
crash. As a side result, a new method is proposed for obtaining regular initial data for Black Hole 
spacetimes, even inside the horizons. 
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I. INTRODUCTION 

The choice of a coordinate system is mandatory when 
finding solutions of Einstein's field equations. In the 3+1 
approach this choice can be made in two steps: 

• The choice of a time coordinate. In geometrical 
terms, this amounts to select a specific foliation of 
spacetime by a family of spacelike hypersurfaces 
(the constant time slices). 

• The choice of the space coordinates. In geometrical 
terms, this amounts to select a specific congruence 
of lines threading every slice of the time foliation 
(the time lines). 

Most of the theoretical developments in the 3+1 ap- 
proach are space- invariant, that is independent of the 
choice of the space coordinates. In this case, one can 
safely select as time lines precisely the normal lines to 
the time slices (normal coordinates), as we will do in 
what follows. In local adapted coordinates, we will write 
then the line element as 

ds 2 = -a 2 dt 2 + jij dx l dx j . (1) 

The choice of a time coordinate, however, is not so sim- 
ple. In numerical simulations one must care about phys- 
ical singularities that could arise dynamically in gravita- 
tional collapse scenarios. Care must be taken to avoid 
the time slices to hit the singularity in a finite amount of 
coordinate time: otherwise the numerical code will crash 
due to the singular terms. 

This is why singularity avoidant slicing conditions are 
widely used in Numerical Relativity. A first example is 
provided by the maximal slicing condition, 



where 

Kij =-^ d t Hi ( 3 ) 

is the extrinsic curvature tensor of the time slices. Maxi- 
mal slicing has been successfully used in Black Hole sim- 
ulations since the early times 0. The behavior of maxi- 
mal slicing in simple cases has been studied for more than 
thirty years Q and is still being studied today 4]. The 
problem with maximal slicing is that it leads to a second 
order elliptic equation for the lapse function a. This can 
be a complication in the current 3D simulations, where 
the need of accuracy is consuming all the available com- 
putational resources. But this is also a complication from 
the theoretical point of view, because the resulting evolu- 
tion system (field equations plus coordinate conditions) 
is then of a mixed elliptic-hyperbolic type. 

A simpler alternative is provided by the generalized 
harmonic condition Q 

d t a = -fa 2 ixK, (4) 

where / can be an arbitrary function of a (/ = 1 in 
the original harmonic slicing 0, Q ) • It provides a first 
order evolution equation for the lapse function which, 
when combined with the field equations, leads to a hyper- 
bolic evolution system (weakly hyperbolic or strongly hy- 
perbolic, depending on the formulation) for non-negative 
choices of /. 

The problem with condition (0J is that it can lead to 
gauge-related pathologies. In Ref. y|, this has been seen 
to occur in Black Hole simulations for some particular 
cases. In Ref. 0, the generic case is considered; starting 
from an analysis of the characteristic speeds, the conclu- 
sion is that there can be two classes of shocks: 

• one which affects just the gauge degrees of freedom. 
A cure is suggested, consisting in the choice 



tr.fi: = 0, 



(2) 



/ = 1 + k/a 2 



(5) 
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{k being just an integration constant). 

• another one which affects even the metric degrees of 
freedom (with c as characteristic speed), for which 
there is no cure. 

The claim concerning the existence of gauge shocks has 
been recently strengthened 01 by providing a purely 
kinematical derivation of the condition JSJ, indepen- 
dently of the field equations. 

These are by no means rhetorical questions. We are 
aware that Einstein's equations can be used for describ- 
ing the propagation of discontinuities arising either from 
the initial or the boundary data. We are talking here 
about something stronger: shocks that could arise dy- 
namically even from smooth initial and boundary data. 
This is a well known phenomenon in the Fluid Dynamics 
domain: the non-linearity of the equations leads in some 
circumstances to the crossing of the characteristic lines 
so that a genuine shock (not just a contact discontinuity) 
appears. Can this really happen in General Relativity?. 
Or, at least, can this happen when using Q as slicing 
condition?. We think that this is a good opportunity to 
address these questions and provide, as far as possible, 
definite answers. 



II. ARE THERE METRIC SHOCKS? 

No, there are not metric shocks. To see why, let us 
start with the space-invariant expression of the corre- 
sponding characteristic speed (light speed) 



Geodesic slicing 




Length 



FIG. 1: The diagram shows a set of light rays as straight lines 
with unit slope in a time- length diagram (geodesic slicing). 
Two hypersurfaces (time slices) of an alternative slicing are 
also shown. These hypersurfaces must be spacelike, so that 
the slope of the tangent vectors is limited by that of the light 
rays. Notice that the space distance among light rays, as 
measured on the time slices, is no longer constant. But notice 
also that the order of the sequence of light rays on each time 
slice is unaltered. 



remains, but condition (J5J) for the proposed cure is to be 
replaced by 



/ = k/a 2 



(7) 



which no longer contains the original harmonic slicing 
(/ = 1) as a particular case. 



dl 



dx % dxi 



dt 



= ± a 



(6) 



(length variation per unit coordinate time). We can al- 
ways consider a geodesic slicing, so that a — 1 and co- 
ordinate time coincides with proper time. In this case, 
as depicted in Fig. ^ the characteristic lines (light rays) 
along any given space direction can not even approach 
one another. 

The question is whether these light rays can cross each 
other in some alternative slicing. As shown in Fig. Q this 
can never happen as long as the time slices are spacelike 
hypersurfaces. Without characteristic crossing, there are 
not genuine shocks (although contact discontinuities can 
appear from initial or boundary data, as discussed in 
the Introduction). Einstein's field equations can then be 
said to be, in the Applied Mathematics jargon, linearly 
degenerate. The case of the gauge degrees of freedom, 
which are not prescribed by the field equations, will be 
discussed in the following sections. 

But let us look back before at Ref. Q, where a 
coordinate-dependent expression for light speed was used 
as the starting point. If we repeat the arguments of 
Ref. 0, but using instead the space-invariant expression 
© as the starting point, we see that the prediction of 
metric shocks disappears. The prediction of gauge shocks 



III. IS THERE A KINEMATICAL REASON FOR 
GAUGE SHOCKS? 

No, there is no kinematical reason for gauge shocks. 
The argument in Ref. [10] starts from the time-slicing 
equation 



1. 



d»(j> d v 4> 



f {~gp"d p <i>d 



Vud v 6 = , 



(8) 



which reduces to the generalized harmonic slicing condi- 
tion (@J in local adapted coordinates, where 



= t. 



(9) 



The argument proceeds then by transforming the quasi- 
linear second-order equation JSJ into a first order system 
and studying the characteristic speeds, without recourse 
to the field equations [To|. 

We will take here a more direct approach, looking for 
the characteristic hypersurfaces of the original equation 
©. To simplify the analysis, we will adapt the general 
coordinates in (JSJ) to the (generic) initial data. Starting 
from a given spacelike slice S, we will choose then as 
initial conditions 



tn 



d k cj) 



0. 



dt</> 



a . 



(10) 
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Notice that this is by no means a restriction on the ini- 
tial data: it is rather the result of adapting the arbi- 
trary spacetime coordinates to any specific choice of ini- 
tial data. Moreover, this kind of adapted coordinates is 
precisely the one that is currently used in numerical sim- 
ulations, where the equation (JHJ is always used in the 
form Q. 

It is trivial now to look for the characteristic hyper- 
surfaces. All second partial derivatives on E can be ob- 
tained from (|10|l with the only exception of the second 
time derivative, which can be obtained from JSJl unless 



I") 



which would then imply that E is a characteristic hyper- 
surface. For any finite choice of /, condition which 
can be written in an invariant way as 



g^d^dvt | B =0, 



(12) 



amounts to the statement that E is a null hypersurface. 
This means that the characteristic speed coincides again 
with light speed. It follows that the discussion of the 
preceding section applies here again to show that there 
can not be any crossing of characteristics and hence no 
genuine shocks can develop. 

To put it in a different way: equation (JSJ) can always 
be used for constructing a new foliation in a given (reg- 
ular) spacetime provided that the resulting slices remain 
spacelike. Then, the time slicing itself will be well de- 
fined, although the behavior of the gauge-related dynam- 
ical fields may depend on the field equations, as we will 
see below. 



IV. ARE THERE DYNAMICAL GAUGE 
SHOCKS? 

There are certainly gauge pathologies of a dynamical 
origin, but not genuine shocks. In order to analyze this 
point, it is convenient to consider the time derivative of 
the generalized harmonic slicing condition Q). This re- 
quires using the field equations to provide the time evo- 
lution of tr K (the arguments in this section will be then 
of a dynamical nature). 

We will choose for simplicity the BSSN formulation [Til 
Il2| . although all the current formalisms share the same 
physical solutions. In the vacuum case, we get after some 
algebra 

(1//) d tt a - A a = - [ a tr (K 2 ) - /' (tr K) 2 ] , (13) 
where we have noted for short 



f = fa 2 



I 



<9< 



(14) 



Equation (|13|l can be interpreted as a generalized wave 
equation for the lapse function. The left-hand-side terms 



correspond to the principal part, and it is clear that the 
characteristic speed is given by 



v G 



±Vf<* 



(15) 



(gauge speed). We derived condition (JJJ by using pre- 
cisely (|15|l as the starting point for the argument pre- 
sented in Ref. Q: it amounts to the requirement of a 
constant gauge speed, that is 



vg 



constant <^> /' = . 



(16) 



The right-hand-side terms in l|13|) can be interpreted 
as non-linear source terms for the evolution of the lapse 
function. Allowing for Q, tr K is proportional to the 
time derivative of a, so that there is a non- linear coupling 
that can lead to runaway (growing without bound, at an 
increasing rate) solutions. 

An extreme example of this gauge pathology, in which 
transport plays not role at all, is provided by space- 
homogeneous and isotropic solutions, like the ones used 
in Cosmology, where all space derivatives vanish and the 
extrinsic curvature is proportional to the metric. One 
has then 

dtt a = /(/' - |) (tr X) 2 = i (/' - |) (d t a) 2 , (17) 

so that a can grow without bound, with a driving force 
proportional to the square of the growing rate, unless 



a 



f = n + 



(n<l). (18) 



Coming back to the general case, it is clear that con- 
dition (|16fl could be interpreted in two alternative ways 

• Ensuring a constant gauge speed, so that no char- 
acteristic crossing can occur. 

• Ensuring the negative sign of the source terms in 
<|13[l . so that no runaway solutions can occur. 

This ambiguity follows from the very nature of the ar- 
gument proposed in Ref. .9j], where the main idea is to 
include the source terms in the discussion about the be- 
havior of characteristic fields. The /-dependent term in 
the right-hand-side of i|13|) would then be considered as a 
non-linear contribution to the transport (left-hand-side) 
terms. 

But we do not see much difference between this partic- 
ular quadratic term and the other ones in the right-hand- 
side of 1)13(1 . independently of their origin. The example 
we presented before shows instead how these other source 
terms actually contribute to the regularity condition l|18|) 
by providing an upper bound which was missing in 116|) . 

Our conclusion is that any pathological behavior of 
the gauge-dependent degrees of freedom, which manifests 
itself as an unbounded growth of either the lapse function 
or its first derivatives, can be interpreted consistently as 
the effect of the non-linear source terms in the evolution 
equations. We will test this conclusion against simple 
numerical Black-Hole simulations in what follows. 
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V. GAUGE PATHOLOGIES IN NUMERICAL 
BLACK-HOLE SIMULATIONS 

We will focus here on three-dimensional numerical sim- 
ulations of a spherically symmetric Black Hole. As far as 
we are interested in the strong field region, the Black Hole 
interior is not excised from the numerical grid. Regular 
initial data are obtained by using a novel 'Free Black 
Hole' approach, which relies on the well known 'no hair' 
theorems (see Appendix I for details). This new approach 
is also being applied independently in the context of the 
BSSN formalism 0. 

We will use in our simulations the first order ver- 
sion of the Z4 evolution formalism [14l Il5| in normal 
coordinates (zero shift). Apart from standard centered 
finite-difference algorithms, we will use at some point the 
MMC method. This is a high-resolution shock-capturing 
(HRSC) method which combines the monotonic-centered 
(MC) flux limiter [3 with Marquina's flux formula [T3 
(see Appendix II for details). The use of these advanced 
HRSC methods is also new in 3D Black Hole simulations, 
although it is crucial for our arguments about shocks. 



A. Slice stretching 

The first stages of the simulation show a collapse of 
the space volume element which, allowing for the singu- 
larity avoidance properties of the gauge conditions l@J, 
translates itself into a collapse of the lapse function a. 
This lapse collapse can be slower or faster, depending on 



FIG. 2: Lapse collapse for the '1+log' slicing : the lapse values 
are shown every Q.5M. The lapse gets very close to zero in 
the innermost region. A standard finite difference algorithm 
is used in a cubic cartesian grid (the space direction in the 
plots corresponds to the grand diagonal of the cube) . Spurious 
oscillations appear at about t — 6M in the zones where the 
slope is changing more abruptly. 
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FIG. 3: Plots of K xx (solid line) and tr K (dotted line) along 
the x axis, corresponding to the 1+log slicing at t = 6M. 
Negative values of K X x correspond to the increasing of -y xx , 
that is a radial stretching of the time slice. Notice that this is 
compatible with an overall collapse pattern, as shown by the 
positivity of tr K: the radial stretching is then compensated 
by the collapse along the angular directions. The profile of 
tr K is smooth, without any sign of shocks. 



our choice of /. We will present here our results for the 
T+log' case 



although no qualitative difference is detected for similar 
choices, like / = 2 for instance. 

Let us see what happens after some time. We see in 
Fig. |21 how high-frequency oscillations start appearing 
in the regions where the lapse slope is changing more 
abruptly. This is a numerical effect: it is the analo- 
gous of the well-known Gibbs phenomenon, which ap- 
pears when the tool one is using (the standard finite 
difference method in our case) is unable to resolve the 
higher frequency modes of some dynamical field. 

In order to see where those annoying high-frequency 
modes came from, let us take a look at behavior of the 
extrinsic curvature, as displayed in Fig. |3| What we see 
there is the well known 'slice stretching' phenomenon: 
the Black Hole is sucking both the numerical grid [Tel ITflj 
and the space slice 0, |2(j, producing as a result steep 
profiles along the radial direction. High-frequency modes 
are the expected to appear in the regions with abrupt 
slope changes in some dynamical field. 

The same kind of spurious oscillations currently ap- 
pears in Computational Fluid Dynamics simulations, in 
regions where genuine shocks are developing. The cure is 
to use advanced HRSC numerical methods that can deal 
with shocks without altering the monotonicity of the dy- 
namical fields. We will use in what follows the MMC 
numerical algorithm, as described in Appendix II, which 
is one of such HRSC methods. We are aware that we 



are not dealing here with real shocks (the tr K profile in 
Fig. 13 is clearly smooth). We are instead dealing with, 
let us say, 'numerical shocks' as far as our numerical grid 
is just unable to resolve the higher frequency modes as- 
sociated with slice stretching. 

The effect of switching to the MMC numerical method 
is definitive, as it can be seen on Fig. 0] No trace of the 
spurious oscillations remains. There is nothing either in 
the lapse or in the tr K profiles that could be interpreted 
as a shock. The conclusion is obvious: slice stretching 
does not produce any real shock, it may produce just 
'numerical shocks' that disappear when improving the 
discretization algorithm. 

B. Runaway solutions 

Let us see now what happens when allowing our sim- 
ulations to proceed for a longer time (about t = 14M in 
our case, but the value will depend on the specific gauge 
choice). We see in Fig. [3] that the spurious oscillations 
do not show up anymore. What we see instead is a sort 
of rebound of the lapse, which is the prelude of a blow 
up that will crash the numerical simulation. We have 
stopped the simulation in Fig. [S] before the crash occurs, 
so that two different interpretations can be tentatively 
considered 

• A gauge shock is forming and starts propagating 
backwards in the region near to p = 2M. 

• The growing of the lapse is triggering a runaway 
solution that will result in a lapse blow up. 



FIG. 4: Same as Figure |5J but using the MMC numerical 
algorithm for the space discretization, as described in Ap- 
pendix II. The spurious oscillations have completely disap- 
peared. Slice stretching is still there (Figure |3] corresponds 
actually to the last plot), but it is no longer a problem for 
numerical simulations. 
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FIG. 5: Same as Figure 0] but with the simulations running 
up to t — 1AM. No spurious oscillations appear anymore. 
However, the lapse collapse stops and a sort of rebound starts 
near the frozen innermost region. The lapse looks just going 
backwards there. The simulation is going to crash by a lapse 
blow up. 



In order to discriminate between these two alterna- 
tives, let us take a look to Fig. where we have plotted 
the same extrinsic curvature components as in Fig. |21 As 
a word of caution, let us remember that the lapse in the 
innermost region is yet collapsed, so that the dynamics 
is frozen there. This means that the features we see in 
the innermost region in Fig. correspond to an earlier 
stage (measured in proper time) than what we see around 
p = 2M, which is the region we are going to analyze now. 

What we see there is just the opposite that what we 
saw in Fig. |3J tr if is now negative, so that we have 
an overall expansion pattern. But the radial component 
(K xx along the x direction) is collapsing. This explains 
why the whole picture in Fig. [S] looks like moving back- 
wards in time in the rebounding region. 

Notice that the tr K profile in Fig. is clearly smooth 
and well-resolved in the region around p = 2M, so that 
the shock interpretation can be excluded. Moreover, 
as far as we are using a HRSC numerical method, no 
blow up would be expected from the presence of shocks. 
In Computational Fluid Dynamics these methods are 
currently used for dealing with genuine shocks without 
crashing because of that. The blow up must be actually 
caused by something much stronger than a mere shock: 
a runaway solution that is triggered by the negative val- 
ues of tr K, which correspond to an overall expansion 
behavior. 

The appearance of runaway solutions is to be expected 
when selecting the T+log' (|19fl or any similar choice 
among the singularity-avoidant family of slicing condi- 
tions Q . If we decompose the extrinsic curvature tensor 
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FIG. 6: Same as Fig. 01 but now for t = 14M, correspond- 
ing to the last plot in Fig. [£] The dynamics in the inner- 
most region is frozen (the lapse is fully collapsed there). In 
the region around p — 2M, however, positive values of K xx 
appear, corresponding to a radial collapse, whereas trK be- 
comes negative, corresponding to a global expansion pattern. 
This is just the behavior opposite to slice stretching, and it 
explains the lapse rebound that can be seen in Fig. [S] Notice 
that the collapse pattern (positive trK values) is recovered 
in the outermost region, where the lapse collapse proceeds as 
expected. 



Kij into its irreducible components 



Aij = 



trK 



K ee trK 



(20) 



(shear tensor and expansion factor, respectively), we see 
that the source terms in the lapse equation 1|13|) can be 
expressed as 



-a (trA 2 



(2-f) K 2 



(21) 



which can easily get a positive sign in collapsing scenarios 
(a small) unless the expansion factor K is kept to low 
values when compared with the shear tensor. 

We conclude that runaway solutions, not gauge shocks, 
can be a real problem in simulations that, like in the 
Black Hole cases, make use of the singularity-avoidant 
slicing conditions (@J. To find a convenient cure for this 
gauge pathology (excising the interior region, using a dy- 
namical shift, selecting other values of /, and so on) goes 
beyond the purpose of this paper. 



Appendix I: Free Black Hole initial data 

Based in the idea that the interior region of a Black 
Hole has no causal physical influence on the exterior one 
(no-hair theorems), we can devise a simple way of ob- 
taining regular initial data for Black Hole spacetimes: 



0.5 



1 . 5 



FIG. 7: Plot of a conformal factor ^ providing 
time-symmetric conformally-flat initial data for a free 
Schwarzschild Black Hole. Values on the x axis correspond 
to the isotropic radial coordinate p, measured in units of M . 
The singular solution in the interior region (dashed line) is 
replaced by a regular one (continuous line). Both expressions 
coincide in the exterior region. Notice that the matching is 
smooth at the apparent horizon (p — Al/2). 



• Solve the energy and momentum constraints for the 
initial data. This currently leads to a singular space 
metric. 

• Replace the metric in the interior regions by any 
smooth extension of the exterior geometry. Energy 
and momentum constraints will be no longer ful- 
filled inside the horizons ('Free Black Hole' data). 

In the simplest case, we can consider time-symmetric 
vacuum initial data with a conformally flat space metric, 
that is 



K. 



ij \t=o 



0. 



lij |t=o 



* 4 <5, 



(LI) 



so that the momentum constraint is trivially satisfied and 
the energy constraint amounts to require that ^ be an 
harmonic function of the Euclidean metric. 

Free initial data can be obtained then for a spherically 
symmetric Black Hole with mass M by starting from the 
singular conformal factor 



* = 1 



M 
V 



(1.2) 



and replacing its value inside the horizon (dotted line in 
Fig.Q by any smooth function (continuous line in Fig. [7| ). 

In the Z4 formalism, energy and momentum constraint 
violations will cause the growing of the supplementary 
fields O and Z^, which non-zero values correspond to a 
departure of true Einstein's solutions. But these fields, 
which are the 3+1 components of the four- vector are 
known to verify |14| 



Rfj.iyZ 1 ' 



0. 



(1.3) 
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so that the spurious non-zero values propagate with light 
speed. The exterior geometry, then, can not be affected 
by constraint violations inside the horizon (at least not 
at the continuum level). 



Appendix II: The MMC high-resolution method 

High-Resolution Shock-Capturing (HRSC) numerical 
methods can be applied to strongly-hyperbolic first-order 
systems of balance laws. The balance-law structure 
means that the evolution equations for the array u of 
dynamical fields can be written as 



d t u + d k F k 



(11.1) 



where both the Flux F k and the Source S terms depend 
algebraically on the fields, but not on their derivatives: 



F fc = F k {u) 



S = S(u) . 



(II.2) 



The system Ijll.lf) will be strongly hyperbolic if and only 
if for every space direction ft, the corresponding charac- 
teristic matrix 



n k A k 



dF k 
du 



(11.3) 



has real eigenvalues v (characteristic speeds) and a com- 
plete set w of eigenvectors |21j . 

The balance law form is specially suited for 

the method-of-lines (MoL) discretization j2^]. In this 



FIG. 8: Starting from the values of the fluxes at the grid 
nodes, a linear reconstruction is done inside every elementary 
cell. Numerical discontinuities appear at every interface (dot- 
ted lines) between the left and right values (arrows and dots, 
respectively) . Notice that the original function was monoton- 
ically decreasing: all the slopes are negative. In the proposed 
reconstruction, however, both the left interface sequence (at 
i — 3/2) and the right interface one (at i + 3/2) show local 
extreme values that break the monotonicity of the original 
function. 



method, there is a clear-cut separation between space 
and time discretization. As a consequence, the source 
terms contribute in a trivial way to the space discretiza- 
tion. The non-trivial contribution comes instead from 
the space derivatives of the Flux terms, which are usu- 
ally discretized as follows 



d k F k 



~Kx [F ^/ 2 



~Kz [F * +1 / 2 



1 *-l/2j 



k-l/2\ 



Ay 1 i+1/2 J-1/2J 

(H.4) 



The half-integer indices correspond to the grid 
' interfaces', which are supposed to be placed halfway be- 
tween neighboring grid nodes. 

Every different way of computing the interface Fluxes 
in terms of the values of the fields at the grid nodes will 
lead to a specific numerical algorithm. A common feature 
of all these algorithms is that the information of every 
grid node is used for providing one-sided predictions at 
the neighbor interfaces in a consistent way. For instance, 
one can define 



i-1/2 



= F, 



Ax . 
— A, 
2 



Ft 



-1/2 



= Fi 



At 

^A,. (II.5) 



Any numerical algorithm must then provide two basic 
elements: 

• A prescription for computing the slopes Aj which 
must be used in the 'reconstruction' process (|II.5|I . 

• A 'Flux formula' that provides a unique value of 
the interface Flux F i+1 / 2 from the one-sided pre- 
dictions {F* 1/2 ,FL +1/2 }. 



FIG. 9: Same as Figure [S] but using the monotonic-centered 
reconstruction. Notice that the interface values are bounded 
now between the neighbor nodes, so that monotonicity is pre- 
served both for the left values (arrows) and for the right ones 
(dots) at every interface (dotted lines). 
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We can see in Fig. [5] an example of the reconstruction 
process for the centered choice 

Af = i(F i+1 - J F i _ 1 ). (II.6) 

The outcome consists in two sequences of one-sided pre- 
dictions at every interface (dots and arrows). Notice that 
the monotonicity of the original Flux sequence is not pre- 
served by the one-sided predictions in the regions with 
abrupt slope changes. This can lead to spurious oscilla- 
tions in the numerical solution, even for smooth dynami- 
cal fields, like the one displayed in Fig. El These spurious 
oscillations are visible in the Black Hole simulation shown 
in Fig.IU 

The MMC method uses instead a monotonicity- 
preserving prescription, which can be obtained by using 
the 'monotonic centered' (MC) slope [l6| 

Af c = minmod{ 2{F l - iVi), 2(F i+1 - F<), A? } 

(H.7) 

instead of (|II.6() . The minmod function is defined as 
usual: 

• If all the arguments have the same sign, then it 
selects the one with smaller absolute value. 

• If one of the arguments has different sign than the 
others, then it is zero. 

In this way, the slopes are limited in order to avoid spu- 
rious oscillations. The rule is that interface values must 
lie between their neighbor node values (see Figure El- 

The second ingredient in the MMC method is the use 
of Marquina's Flux formula [17J . The prescription to cal- 
culate the interface values F i+1 n can be stated as follows 
(we have simplified the original formula by adapting it to 
the quasilinear evolution systems one gets from Einstein's 
field equations): 



• Decompose both one-sided predictions (F L , F R ) 
as linear combinations of the set of characteristic 
fields w. Note that the coefficients in these combi- 
nations are not necessarily constant: we must use 
in general a different set of values on each side of 
the interface. 

• Project the forward prediction F L , by suppressing 
any components w corresponding to negative char- 
acteristic speeds (forward projection Fp). 

• Project the backward prediction F R , by suppress- 
ing any components w corresponding to positive 
characteristic speeds (backward projection Fb). 

• Add these upwind-projected values at the given in- 
terface, that is 

F m/2 = Fjr + Fb ■ {U.8) 

One is taking in this way the positive speed components 
from the previous cell and the negative speed components 
from the next one. Notice that we must assume for con- 
sistency a definite sign for the characteristic speed (just 
the sign, not even the value). If the sign changes between 
both sides of the given interface (when using a super- 
luminal shift, for instance), then another combination 
must be taken instead of the simple upwind projections 
presented here (see Ref. 0] for the details). 
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